#Democratic Reputations in Crisis and War: replication data
#November 20, 2021
#Josh Kertzer

######### Block A Replication II 112021.R: this R file contains the code necessary to replicate the analysis in the main text - run Block A Replication I 112021.R first!
# All of the following analyses were carried out using R version 4.1.2 GUI 1.77 (Big Sur ARM build 8007) on an M1 Pro Macbook Pro running MacOS Monterey. 

#### Table 2: Mass Public Samples

#Column for Israel 1
c1 <- c(mean(isr1$Male, na.rm=TRUE), table(isr1$ageCat)/sum(table(isr1$ageCat)), table(isr1$educCat4)/sum(table(isr1$educCat4)), mean(isr1$ideo1, na.rm=TRUE), mean(isr1$milAssert1, na.rm=TRUE), mean(isr1$intTrust1, na.rm=TRUE), mean(isr1$hawk, na.rm=TRUE), NA, NA, nrow(isr1))

#Column for Israel 2
c2 <- c(mean(isr2$Male, na.rm=TRUE), table(isr2$ageCat)/sum(table(isr2$ageCat)), table(isr2$educCat4)/sum(table(isr2$educCat4)), mean(isr2$ideo1, na.rm=TRUE), mean(isr2$milAssert1, na.rm=TRUE), NA, mean(isr1$hawk, na.rm=TRUE), NA, NA, nrow(isr2))

#Column for ROK
c3 <- c(mean(rok$male, na.rm=TRUE), table(rok$ageCat)/sum(table(rok$ageCat)), table(rok$educCat4)/sum(table(rok$educCat4)), mean(rok$ideo1, na.rm=TRUE), NA, NA, NA, mean(rok$interestGov1, na.rm=TRUE), mean(rok$interestFP1, na.rm=TRUE), nrow(rok))

#Column for UK 1
c4 <- c(mean(uk1$male, na.rm=TRUE), table(uk1$ageCat)/sum(table(uk1$ageCat)), table(uk1$educCat4)/sum(table(uk1$educCat4)), mean(uk1$ideo1, na.rm=TRUE), NA, NA, NA, mean(uk1$interestGov1, na.rm=TRUE), mean(uk1$interestFP1, na.rm=TRUE), nrow(uk1))

#Column for UK 2
c5 <- c(mean(uk2$male, na.rm=TRUE), table(uk2$ageCat)/sum(table(uk2$ageCat)), table(uk2$educCat4)/sum(table(uk2$educCat4)), mean(uk2$ideo1, na.rm=TRUE), NA, NA, NA, mean(uk2$interestGov1, na.rm=TRUE), mean(uk2$interestFP1, na.rm=TRUE), nrow(uk2))

#Column for USA
c6 <- c(mean(usa$Male, na.rm=TRUE), table(usa$ageCat)/sum(table(usa$ageCat)), table(usa$educCat4)/sum(table(usa$educCat4)), mean(usa$ideo1, na.rm=TRUE), mean(usa$milAssert1, na.rm=TRUE), mean(usa$intTrust, na.rm=TRUE), NA, mean(usa$interestGov1, na.rm=TRUE), mean(usa$interestFA1, na.rm=TRUE), nrow(usa))

tab2 <- round(cbind(c1, c2, c3, c4, c5, c6), digits=2) 
colnames(tab2) <- c("Israel I", "Israel II", "Korea", "UK I", "UK II", "USA")
rownames(tab2) <- c("Male", "Age: < 25", "25-34", "35-44", "45-54", ">55", "Education: High School or less", "Some college", "College/university", "Postgraduate", "Ideology", "Mil. Assert", "Int Trust", "Arab-Israeli conflict", "Pol. Interest", "Interest in FP", "N")

stargazer(tab2, digits=2)

#Total N: 8855
sum(tab2[nrow(tab2),])

#### Figure 3: Democratic reputations

#Rescale outcome measures, render them interpretable on a probability scale
#For crises:
isr2$resolve2 <- isr2$DV_A1/100
isr2$publicPacifist2 <- rescale(isr2$publicPacifist)
isr2$elitePacifist2 <- rescale(isr2$elitePacifist)
isr2$casualtySensitive2 <- recode(isr2$casualtySensitive, binarize=TRUE, x=3, x2=4)
isr2$finSensitive2 <- recode(isr2$finSensitive, binarize=TRUE, x=3, x2=4)
isr2$turnoverDiff2 <- recode(isr2$turnoverDiff,binarize=TRUE, x=-4, x2=-1)
isr2$turnoverDiffC2 <- isr2$turnoverDiff2
isr2$turnoverDiffC2[which(isr2$mech!="Closed_Crisis")] <- NA
isr2$pubSupportDiff2 <- recode(isr2$pubSupportDiff,binarize=TRUE, x=0,x2=4)
isr2$pubSupportDiffC2 <- isr2$pubSupportDiff2
isr2$pubSupportDiffC2[which(isr2$mech!="Closed_Crisis")] <- NA
isr2$irregularExit2 <- recode(as.numeric(isr2$llFate),binarize=TRUE, x=3,x2=5)
isr2$irregularExitC2 <- isr2$irregularExit2
isr2$irregularExitC2[which(isr2$mech!="Closed_Crisis")] <- NA
isr2$followThrough2 <- recode(isr2$followThrough, binarize=TRUE, x=4,x2=5)
#For wars:
isr2$win2 <- isr2$DV_A2/100
isr2$selectWin2 <- recode(isr2$selectWin, binarize=TRUE, x=4,x2=5)
isr2$turnoverDiffW2 <- isr2$turnoverDiff2
isr2$turnoverDiffW2[which(isr2$mech!="Closed_War")] <- NA
isr2$pubSupportDiffW2 <- isr2$pubSupportDiff2
isr2$pubSupportDiffW2[which(isr2$mech!="Closed_War")] <- NA
isr2$irregularExitW2 <- isr2$irregularExit2
isr2$irregularExitW2[which(isr2$mech!="Closed_War")] <- NA
isr2$reliableAllies2 <- recode(isr2$reliableAllies, binarize=TRUE, x=4,x2=5)
isr2$training2 <- recode(isr2$training, binarize=TRUE, x=4,x2=5)
isr2$morale2 <- recode(isr2$morale, binarize=TRUE, x=4, x2=5)

## Parallel processing (to speed up the analysis)
cl <- makeCluster(20,"SOCK") #Modify for the number of cores available on your computer
#For bootstrapping
B <- 1500
set.seed(43215)

bootA <- function(...){
	j <- sample(1:nrow(isr2), replace=TRUE)
	temp.out <- rep(NA, length(x.vec))
	X <- isr2$A[j]
	for (k in 1:length(x.vec)){
		Y <- isr2[j,x.vec[k]]
		temp.out[k] <- mean(Y[X==1],na.rm=TRUE) - mean(Y[X==0],na.rm=TRUE)
	}
	return(temp.out)
}

clusterSetRNGStream(cl, 43215) #Set seed for parallel processing

## Panel a: Democratic reputations in crises

x.vec <- c("resolve2","publicPacifist2","elitePacifist2","casualtySensitive2","finSensitive2", "turnoverDiffC2", "pubSupportDiffC2", "irregularExitC2", "followThrough2") #Vector of variable names

clusterExport(cl, list("isr2", "x.vec", "bootA"), envir=environment())
mat.c <- parSapply(cl,rep(1,B),bootA) #Obtain distributions

mat.c2 <- t(mat.c)
c.vec.lab <- c("Stand firm in crises", "Dovish public", "Dovish leaders", "Casualty sensitive", "Financially sensitive", "Leader turnover", "Maintain public support", "Irregular exit", "Follow through") #Crisis plot labels

colnames(mat.c2) <- c.vec.lab

df.c <- gather(as.data.frame(mat.c2))
df.c$key <- factor(df.c$key, levels=rev(c.vec.lab))

dev.new(height=5, width=5)
ggplot(df.c, aes(y=key)) + geom_density_ridges(aes(x=value)) + theme_bw() + xlab("Effect of democracy in crises") + ylab("") + geom_vline(xintercept=0, lty=3) + xlim(-0.5,0.5)
#Add additional labels to plot in post-production

## Panel b: Democratic reputations in wars

x.vec <- c("win2", "selectWin2", "turnoverDiffW2", "pubSupportDiffW2", "irregularExitW2", "reliableAllies2", "training2", "morale2")

clusterExport(cl, list("isr2", "x.vec", "bootA"), envir=environment())
mat.w <- parSapply(cl,rep(1,B),bootA)

mat.w2 <- t(mat.w)
w.vec.lab <- c("Win in war", "Pr(Victory|Initiate)", "Leader turnover",  "Maintain public support", "Irregular exit", "Allies", "Training", "Morale") #War plot labels
colnames(mat.w2) <- w.vec.lab

df.w <- gather(as.data.frame(mat.w2))
df.w$key <- factor(df.w$key, levels=rev(w.vec.lab))

dev.new(height=5, width=5)
ggplot(df.w, aes(y=key)) + geom_density_ridges(aes(x=value)) + theme_bw() + xlab("Effect of democracy in war") + ylab("") + geom_vline(xintercept=0, lty=3) + xlim(-0.6,0.4)
#Add additional labels to plot in post-production

## Descriptive statistics in text

round(apply(mat.w2,2,mean),digits=2)[1] #Democracies seen as 8% more likely to win in war
round(apply(mat.c2,2,mean), digits=2) #Democracies seen as 2% less likely to win in crises, 10% more likely to have public with dovish preferences, 16% more likely to have leaders with dovish preferences, 37% more likely to be casualty sensitive, 23% more likely to financially sensitive, 3% more likely to lose power, 2% less likely to maintain public support, 35% less likely to experience irregular exit, 7% less likely to follow through

round(apply(mat.w2,2,mean),digits=2)[c(2,6:8)] #Democracies seen as 16% more likely to win if it initiated (vs non-democracies), 18% more reliable allies, 14% better training, 23% better morale

round(apply(mat.w2,2,mean),digits=2)[c(3:5)] #Democracies seen as 5% more likely to lose power, 6% less likely to maintain public support, 52% less likely to experience irregular exit

#### Figure 4: Average treatment effects and difference-in-differences

## Panel a: Density plot of ATEs

set.seed(43215)
B <- 1500

###Calculate UK public I
boot.uk1.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(uk1), replace=TRUE)
  temp <- uk1[j,]
  boot.uk1.1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

boot.uk1.2 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(uk1), replace=TRUE)
  temp <- uk1[j,]
  boot.uk1.2[i] <- mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE) 
}

### Calculate UK public II

boot.uk2.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(uk2), replace=TRUE)
  temp <- uk2[j,]
  boot.uk2.1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

boot.uk2.2 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(uk2), replace=TRUE)
  temp <- uk2[j,]
  boot.uk2.2[i] <- mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE) 
}


### Calculate ROK public
boot.rok.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(rok), replace=TRUE)
  temp <- rok[j,]
  boot.rok.1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

boot.rok.2 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(rok), replace=TRUE)
  temp <- rok[j,]
  boot.rok.2[i] <- mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE) 
}

#Calculate Knesset
knes.dat1 <- na.omit(as.data.frame(cbind(DV_A1=knes$DV_A1, A=knes$A))) #Drop missingness
boot.k1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(knes), replace=TRUE)
  temp <- knes.dat1[j,]
  boot.k1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

knes.dat2 <- na.omit(as.data.frame(cbind(DV_A2=knes$DV_A2, A=knes$A))) #Drop missingness
boot.k2 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(knes.dat2), replace=TRUE)
  temp <- knes.dat2[j,]
  boot.k2[i] <- mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE) 
}

#Calculate Israeli Public I
boot2.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(isr1), replace=TRUE)
  temp <- isr1[j,]
  boot2.1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

boot2.2 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(isr1), replace=TRUE)
  temp <- isr1[j,]
  boot2.2[i] <- mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE) 
}


### Calculate Israeli Public II
boot.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(isr2), replace=TRUE)
  temp <- isr2[j,]
  boot.1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

boot.2 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(isr2), replace=TRUE)
  temp <- isr2[j,]
  boot.2[i] <- mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE) 
}

### Calculate US
usa2 <- usa[which(usa$RegimeB=="Dictatorship"),] #Focus on comparable subset
boot.usa1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(usa2), replace=TRUE)
  temp <- usa2[j,]
  boot.usa1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

boot.usa2 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(usa2), replace=TRUE)
  temp <- usa2[j,]
  boot.usa2[i] <- mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE) 
}


x_2rev <- data.frame(resolve=c(boot.uk1.1, boot.uk2.1, boot.rok.1, boot.k1, boot2.1, boot.1, boot.usa1), effective=c(boot.uk1.2, boot.uk2.2, boot.rok.2, boot.k2, boot2.2, boot.2, boot.usa2), z=rep(c("United Kingdom I", "United Kingdom II", "Republic of Korea", "Israel: Knesset", "Israel I", "Israel II", "USA"), each=B))

dev.new(height=5,width=4.5)
ggplot(x_2rev) + geom_density(aes(x=resolve, y=-..density.., fill=z), alpha=0.55) + geom_density(aes(x=effective, y=..density..,fill=z), alpha=0.55) + xlim(-20,25) + theme_bw() + scale_y_continuous(breaks=c(-.2,.2), labels=c("...in crises", "...in wars")) + geom_vline(xintercept=0, linetype=3, alpha=0.5, show.legend=FALSE) + labs(x="Effect of democracy", y="") + scale_fill_brewer(palette="Set1") + guides(fill="none") + coord_flip()

### Panel b: Difference-in-difference
set.seed(43215)
B <- 1500

###Calculate UK public I
diff.uk1.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(uk1), replace=TRUE)
  temp <- uk1[j,]
  diff.uk1.1[i] <- (mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE)) - (mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE)) 
}

### Calculate UK public II

diff.uk2.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(uk2), replace=TRUE)
  temp <- uk2[j,]
  diff.uk2.1[i] <- (mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE)) - (mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE)) 
}

### Calculate ROK public
diff.rok.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(rok), replace=TRUE)
  temp <- rok[j,]
  diff.rok.1[i] <- (mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE)) - (mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE)) 
}


#Calculate Knesset
diff.k1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(knes), replace=TRUE)
  temp <- knes[j,]
  diff.k1[i] <- (mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE)) - (mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE)) 
}


#Calculate Israeli Public I
boot2.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(isr1), replace=TRUE)
  temp <- isr1[j,]
  boot2.1[i] <- (mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE)) - (mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE))
}

### Calculate Israeli Public II
diff.1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(isr2), replace=TRUE)
  temp <- isr2[j,]
  diff.1[i] <- (mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE)) - (mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE)) 
}

### Calculate US MTurk

diff.usa1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(usa2), replace=TRUE)
  temp <- usa2[j,]
  diff.usa1[i] <- (mean(temp$DV_A2[temp$A==1], na.rm=TRUE) - mean(temp$DV_A2[temp$A==0], na.rm=TRUE)) - (mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE)) 
}

x_3 <- data.frame(y=c(diff.uk1.1, diff.uk2.1, diff.rok.1, diff.k1, boot2.1, diff.1, diff.usa1), z=rep(c("United Kingdom I","United Kingdom II", "Republic of Korea","Israel: Knesset", "Israel I", "Israel II", "USA"), each=B))

dev.new(height=5,width=6.2)
ggplot(x_3) + geom_density(aes(x=y, y=..density.., fill=z), alpha=0.55) + scale_fill_brewer(guide=guide_legend(title="Sample"), palette="Set1") + labs(x="Difference-in-difference: effect of democracy in war - crisis", y="") + scale_y_continuous(breaks=NULL) + theme_bw() + geom_vline(xintercept=0, linetype=3, alpha=0.5, show.legend=FALSE) #+ coord_flip()

## Descriptive statistics mentioned in text

library(dplyr)
x_2rev %>% group_by(z) %>% summarize(crises=mean(resolve), war=mean(effective))
#Two samples where democracies are seen as being at a disadvantage
#Knesset (-5.9%) and Israel I (-4.6%)
#Democratic advantage in war particularly pronounced in Knesset sample: 14.4%

## Footnote 14:

round(summary(lm(DV_A1 ~ A*currentMK, data=knes))$coef[4,4],digits=2) #p < 0.13 for interaction with current MK, resolve
round(summary(lm(DV_A1 ~ A*milAssert1, data=knes))$coef[4,4],digits=2) #p < 0.40 for interaction with current MK, resolve
round(summary(lm(DV_A2 ~ A*currentMK, data=knes))$coef[4,4],digits=2) #p < 0.71 for interaction with current MK, military effectiveness
round(summary(lm(DV_A2 ~ A*milAssert1, data=knes))$coef[4,4],digits=2) #p < 0.41 for interaction with current MK, military effectiveness